1 Goal of the script

This script plots all variables to see which ones should be used for further analysis.

dir_in <- "analysis_before.after/derived_data/"
dir_out <- "analysis_before.after/plots"

Raw data must be located in ~/analysis_before.after/derived_data/.
Formatted data will be saved in ~/analysis_before.after/plots.

1.1 The knit directory for this script is the project directory.

2 Load packages

pack_to_load <- c("R.utils", "tools", "ggplot2", "doBy", "tidyverse", "patchwork", 
                "ggsci", "ggfortify")
sapply(pack_to_load, library, character.only = TRUE, logical.return = TRUE)
Warning: package 'R.utils' was built under R version 4.1.3
Warning: package 'doBy' was built under R version 4.1.3
Warning: package 'tibble' was built under R version 4.1.3
Warning: package 'tidyr' was built under R version 4.1.3
Warning: package 'readr' was built under R version 4.1.3
Warning: package 'dplyr' was built under R version 4.1.3
Warning: package 'ggfortify' was built under R version 4.1.3
  R.utils     tools   ggplot2      doBy tidyverse patchwork     ggsci ggfortify 
     TRUE      TRUE      TRUE      TRUE      TRUE      TRUE      TRUE      TRUE 

3 Get name, path and information of the file

data_file <- list.files(dir_in, pattern = "\\.Rbin$", full.names = TRUE)
md5_in <- md5sum(data_file)
info_in <- data.frame(file = basename(names(md5_in)), checksum = md5_in, 
           row.names = NULL)

4 Load data into R object

imp_data <- loadObject(data_file)
str(imp_data)
'data.frame':   144 obs. of  43 variables:
 $ Sample                  : chr  "FLT4-10" "FLT4-10" "FLT4-10" "FLT4-10" ...
 $ Spot                    : chr  "A" "A" "B" "B" ...
 $ Cycle                   : Factor w/ 2 levels "before","2000": 2 1 2 1 2 1 2 1 2 1 ...
 $ Raw.material            : Factor w/ 2 levels "flint","lydite": 1 1 1 1 1 1 1 1 1 1 ...
 $ Contact.material        : chr  "bone plate" "bone plate" "bone plate" "bone plate" ...
 $ Acquisition.date.time   : chr  "08.04.2021 13:26" "08.05.2021 07:37" "08.04.2021 14:07" "08.05.2021 07:56" ...
 $ Analysis.date           : chr  "04.10.2022" "04.10.2022" "04.10.2022" "04.10.2022" ...
 $ Analysis.time           : 'times' num  0.635 0.636 0.636 0.636 0.637 ...
  ..- attr(*, "format")= chr "h:m:s"
 $ Sq                      : num  2681 1772 358 425 1219 ...
 $ Ssk                     : num  -0.0173 1.8608 3.0106 2.1596 0.7077 ...
 $ Sku                     : num  1.99 13.5 24.97 11.13 3.21 ...
 $ Sp                      : num  7727 11591 3313 2962 4083 ...
 $ Sv                      : num  5870 5604 1272 1134 2779 ...
 $ Sz                      : num  13597 17195 4585 4096 6862 ...
 $ Sa                      : num  2328 1083 224 271 938 ...
 $ Smr                     : num  0.389 0.277 0.45 0.6 1.56 ...
 $ Smc                     : num  3277 1301 277 408 1865 ...
 $ Sxp                     : num  4666 3448 559 568 1783 ...
 $ Sal                     : num  24.93 16.8 8.93 8.43 18.27 ...
 $ Str                     : num  NA 0.445 0.33 0.255 NA ...
 $ Std                     : num  176 130 168 150 176 ...
 $ Sdq                     : num  0.287 0.329 0.123 0.195 0.234 ...
 $ Sdr                     : num  3.751 4.435 0.735 1.559 2.503 ...
 $ Vm                      : num  0.0783 0.2056 0.0371 0.0571 0.0662 ...
 $ Vv                      : num  3.355 1.507 0.314 0.466 1.931 ...
 $ Vmp                     : num  0.0783 0.2056 0.0371 0.0571 0.0662 ...
 $ Vmc                     : num  2.891 1.039 0.231 0.225 1.052 ...
 $ Vvc                     : num  3.172 1.289 0.281 0.431 1.839 ...
 $ Vvv                     : num  0.1833 0.2175 0.0336 0.0345 0.093 ...
 $ Maximum.depth.of.furrows: num  5229 7056 1817 1907 4086 ...
 $ Mean.depth.of.furrows   : num  1236 1178 355 749 1546 ...
 $ Mean.density.of.furrows : num  3460 3531 4040 3956 3547 ...
 $ Isotropy                : num  NA 44.5 33 25.5 NA ...
 $ First.direction         : num  1.80e+02 1.36e+02 1.80e+02 3.89e-03 1.80e+02 ...
 $ Second.direction        : num  135.7181 0.0166 137.2588 89.9832 153.6998 ...
 $ Third.direction         : num  90 130 90 137 90 ...
 $ Texture.isotropy        : num  NA 42.5 48.5 25.7 NA ...
 $ epLsar                  : num  0.00572 0.00297 0.00249 0.00513 0.00482 ...
 $ NewEplsar               : num  0.0202 0.0186 0.0187 0.0196 0.0196 ...
 $ Asfc                    : num  4.47 6.44 1.44 3.02 3.62 ...
 $ Smfc                    : num  7.64 4.07 2.62 3.58 7.84 ...
 $ HAsfc9                  : num  0.313 0.556 0.257 0.766 0.191 ...
 $ HAsfc81                 : num  0.584 0.846 0.523 1.979 0.604 ...
 - attr(*, "comment")= Named chr [1:42] "µm" "points" "µm" "µm" ...
  ..- attr(*, "names")= chr [1:42] "Axis length - X" "Axis size - X" "Axis spacing - X" "Axis length - Y" ...

The imported file is: “~/analysis_before.after/derived_data/before.after.Rbin”


5 Prepare variables

5.1 Define numeric variables

num.var <- 9:length(imp_data)

The following variables will be used:

[9] Sq
[10] Ssk
[11] Sku
[12] Sp
[13] Sv
[14] Sz
[15] Sa
[16] Smr
[17] Smc
[18] Sxp
[19] Sal
[20] Str
[21] Std
[22] Sdq
[23] Sdr
[24] Vm
[25] Vv
[26] Vmp
[27] Vmc
[28] Vvc
[29] Vvv
[30] Maximum.depth.of.furrows
[31] Mean.depth.of.furrows
[32] Mean.density.of.furrows
[33] Isotropy
[34] First.direction
[35] Second.direction
[36] Third.direction
[37] Texture.isotropy
[38] epLsar
[39] NewEplsar
[40] Asfc
[41] Smfc
[42] HAsfc9
[43] HAsfc81

6 Calculate and plot the absolute difference for each sample and each parameter

6.1 Sorted by contact material

6.1.1 pig skin & skin pad

# add another coloum which concatenates Sample and Spot 
imp_data$ID <- paste(imp_data$Sample, imp_data$Spot)

diff_data <- imp_data[-(6:8)] %>% 
  pivot_wider(names_from = Cycle, values_from = all_of(names(imp_data)[num.var]))


output_list <- vector(mode = "list", length = length(names(imp_data)[num.var]))
names(output_list) <- paste0(names(imp_data)[num.var], ".diff")

for (i in seq_along(num.var)){
  temp <- select(diff_data, starts_with(paste0(names(imp_data[num.var])[i], "_")))
  output_list[[i]] <- unlist(select(temp, contains("2000")) - select(temp, 
                      contains("before")))
}

var.diff <- as.data.frame(output_list)

# create a data frame with the differences from before and after as values 
difference <- as.data.frame(cbind(select(diff_data, Sample:ID), var.diff))
row.names(difference) <- NULL


# define numeric variable 
num.var2 <- 6:length(difference)
# the following variables will be used: 
for (j in num.var2) cat("[",j,"] ", names(difference)[j], "\n", sep="")
[6] Sq.diff
[7] Ssk.diff
[8] Sku.diff
[9] Sp.diff
[10] Sv.diff
[11] Sz.diff
[12] Sa.diff
[13] Smr.diff
[14] Smc.diff
[15] Sxp.diff
[16] Sal.diff
[17] Str.diff
[18] Std.diff
[19] Sdq.diff
[20] Sdr.diff
[21] Vm.diff
[22] Vv.diff
[23] Vmp.diff
[24] Vmc.diff
[25] Vvc.diff
[26] Vvv.diff
[27] Maximum.depth.of.furrows.diff
[28] Mean.depth.of.furrows.diff
[29] Mean.density.of.furrows.diff
[30] Isotropy.diff
[31] First.direction.diff
[32] Second.direction.diff
[33] Third.direction.diff
[34] Texture.isotropy.diff
[35] epLsar.diff
[36] NewEplsar.diff
[37] Asfc.diff
[38] Smfc.diff
[39] HAsfc9.diff
[40] HAsfc81.diff
for (j in num.var2){
  
  # get the min/max range of the data set
  range_var <- range(difference[j])  
  # plot pig skin 
  p_dif_skin <- ggplot(data = difference[grep("pig skin", difference[["Contact.material"]]),
                ], aes_string(x = "Spot", names(difference)[j], colour = "Sample")) + 
                  geom_jitter(size = 3, width = 0.25) +
                coord_cartesian(ylim = range_var) +
                  facet_wrap(Sample ~ Contact.material, nrow = 2) + 
                labs(title = "natural contact material") + 
                ylab(names(difference)[j]) + xlab(NULL) +
                  labs(y = gsub("\\.", " ", names(difference)[j])) +
                  scale_colour_futurama() +
                  theme_classic()+
                theme(legend.position = "none") 

   # plot skin pad 
   p_dif_skin.pad <- ggplot(data = difference[grep("skin pad",
                     difference[["Contact.material"]]), ], aes_string(x = "Spot",    
                     names(difference)[j], 
                       colour = "Sample")) + 
                       geom_jitter(size = 3, width = 0.25) +
                     
                     coord_cartesian(ylim = range_var) +
                       facet_wrap(Sample ~ Contact.material, nrow = 2) + 
                     labs(title = "artificial contact material") + 
                     ylab(names(difference)[j]) + xlab(NULL) +
                       labs(y = gsub("\\.", " ", names(difference)[j])) +
                       scale_colour_futurama() +
                       theme_classic() +
                     theme(legend.position = "none") 

    # combine both plots 
  p_dif.skin <- p_dif_skin + p_dif_skin.pad + plot_layout(width = c(3/6, 3/6))  
  print(p_dif.skin)

  # save to PDF
    file_out <- paste0(file_path_sans_ext(info_in[["file"]]), "_skin_", names(difference)[j],
                ".pdf")
    ggsave(filename = file_out, plot = p_dif.skin, path = dir_out, device = "pdf", 
           width = 300, height = 170, units = "mm")

}

Warning: Removed 5 rows containing missing values (geom_point).
Removed 5 rows containing missing values (geom_point).
Removed 5 rows containing missing values (geom_point).
Removed 5 rows containing missing values (geom_point).

Warning: Removed 5 rows containing missing values (geom_point).
Removed 5 rows containing missing values (geom_point).
Removed 5 rows containing missing values (geom_point).
Removed 5 rows containing missing values (geom_point).

Warning: Removed 1 rows containing missing values (geom_point).
Warning: Removed 1 rows containing missing values (geom_point).

7 cow scapula & bone plate

for (j in num.var2){
 # get the min/max range of the data set
 range_var <- range(difference[j])  
  
  # plot cow scapula 
  p_dif_bone <- ggplot(data = difference[grep("cow scapula",   
                difference[["Contact.material"]]), ], aes_string(x = "Spot", 
                names(difference)[j], colour = "Sample")) + 
                  geom_jitter(size = 3, width = 0.25) +
                coord_cartesian(ylim = range_var) +
                  facet_wrap(Sample ~ Contact.material, nrow = 2) + 
                labs(title = "natural contact material") + 
                ylab(names(difference)[j]) + xlab(NULL) +
                  labs(y = gsub("\\.", " ", names(difference)[j])) +
                  scale_colour_futurama() +
                  theme_classic()+
                theme(legend.position = "none") 
 
   # plot bone plate 
   p_dif_bone.plate <- ggplot(data = difference[grep("bone plate", 
                       difference[["Contact.material"]]), ], aes_string(x = "Spot", 
                       names(difference)[j], colour = "Sample")) + 
                         geom_jitter(size = 3, width = 0.25) +
                       coord_cartesian(ylim = range_var) +
                         facet_wrap(Sample ~ Contact.material, nrow = 2) + 
                       labs(title = "artificial contact material") + 
                       ylab(names(difference)[j]) + xlab(NULL) +
                         labs(y = gsub("\\.", " ", names(difference)[j])) +
                         scale_colour_futurama() +
                         theme_classic()+
                       theme(legend.position = "none") 
 
    # combine both plots 
  p_dif.bone <- p_dif_bone + p_dif_bone.plate + plot_layout(width = c(3/6, 3/6))  
  print(p_dif.bone)

  # save to PDF
    file_out <- paste0(file_path_sans_ext(info_in[["file"]]), "_bone_", names(difference)[j],
                ".pdf")
    ggsave(filename = file_out, plot = p_dif.bone, path = dir_out, device = "pdf", 
           width = 300, height = 170, units = "mm")

}

Warning: Removed 2 rows containing missing values (geom_point).
Warning: Removed 4 rows containing missing values (geom_point).
Warning: Removed 2 rows containing missing values (geom_point).
Warning: Removed 4 rows containing missing values (geom_point).

Warning: Removed 2 rows containing missing values (geom_point).
Removed 4 rows containing missing values (geom_point).
Warning: Removed 2 rows containing missing values (geom_point).
Warning: Removed 4 rows containing missing values (geom_point).

Warning: Removed 1 rows containing missing values (geom_point).
Warning: Removed 3 rows containing missing values (geom_point).
Warning: Removed 1 rows containing missing values (geom_point).
Warning: Removed 3 rows containing missing values (geom_point).

Warning: Removed 1 rows containing missing values (geom_point).
Warning: Removed 1 rows containing missing values (geom_point).

8 Calculate and plot the absolute difference for each sample and each parameter as mean value

8.1 sorted by contact material

8.1.1 cow scapula & bone plate

# compute the mean of the three spots per sample 
mean.sample <- summaryBy(.~ Sample + Raw.material + Contact.material, 
               data = difference, FUN = mean)
# define new num.var for mean_cycle 
num.var3 <- 4:length(mean.sample)


# the following variables will be used: 
for (k in num.var3) cat("[",k,"] ", names(difference)[k], "\n", sep="")
[4] Contact.material
[5] ID
[6] Sq.diff
[7] Ssk.diff
[8] Sku.diff
[9] Sp.diff
[10] Sv.diff
[11] Sz.diff
[12] Sa.diff
[13] Smr.diff
[14] Smc.diff
[15] Sxp.diff
[16] Sal.diff
[17] Str.diff
[18] Std.diff
[19] Sdq.diff
[20] Sdr.diff
[21] Vm.diff
[22] Vv.diff
[23] Vmp.diff
[24] Vmc.diff
[25] Vvc.diff
[26] Vvv.diff
[27] Maximum.depth.of.furrows.diff
[28] Mean.depth.of.furrows.diff
[29] Mean.density.of.furrows.diff
[30] Isotropy.diff
[31] First.direction.diff
[32] Second.direction.diff
[33] Third.direction.diff
[34] Texture.isotropy.diff
[35] epLsar.diff
[36] NewEplsar.diff
[37] Asfc.diff
[38] Smfc.diff
for (k in num.var3){
  # get the min/max range of the data set
  range_var <- range(mean.sample[k])  
  # plot cow scapula
  bone.mean <- ggplot(data = mean.sample[grep("cow scapula", 
               mean.sample[["Contact.material"]]), ], aes_string(x = "Sample", 
               names(mean.sample)[k], colour = "Sample")) + 
               geom_point(size = 4) +
               coord_cartesian(ylim = range_var) +
                 facet_wrap(Sample ~ Contact.material, nrow = 2) + 
               labs(title = "natural contact material") + 
               ylab(names(mean.sample)[j]) + xlab(NULL) +
                 labs(y = gsub("\\.", " ", names(mean.sample)[k])) +
                 scale_colour_futurama() +
                 theme_classic() +
               theme(legend.position = "none") +
               theme(axis.text.x = element_blank(),axis.ticks = element_blank())  
                
   # plot bone plate 
   plate.mean <- ggplot(data = mean.sample[grep("bone plate",
                 mean.sample[["Contact.material"]]), ], aes_string(x = "Sample",
                 names(mean.sample)[k], colour = "Sample")) + 
                   geom_point(size = 4) +
                 coord_cartesian(ylim = range_var) +
                   facet_wrap(Sample ~ Contact.material, nrow = 2) + 
                 labs(title = "artificial contact material") + 
                 ylab(names(mean.sample)[j]) + xlab(NULL) +
                   labs(y = gsub("\\.", " ", names(mean.sample)[k])) +
                   scale_colour_futurama() +
                   theme_classic() +
                 theme(legend.position = "none") +
                 theme(axis.text.x = element_blank(),axis.ticks = element_blank())  


    # combine both plots 
  bone.plate.mean <- bone.mean + plate.mean + plot_layout(width = c(3/6, 3/6))  
  print(bone.plate.mean)

  # save to PDF
    file_out <- paste0(file_path_sans_ext(info_in[["file"]]), "_bone_", names(mean.sample)[k],
                ".pdf")
    ggsave(filename = file_out, plot = bone.plate.mean, path = dir_out, device = "pdf", 
           width = 300, height = 170, units = "mm")

} 

Warning: Removed 2 rows containing missing values (geom_point).
Warning: Removed 3 rows containing missing values (geom_point).
Warning: Removed 2 rows containing missing values (geom_point).
Warning: Removed 3 rows containing missing values (geom_point).

Warning: Removed 2 rows containing missing values (geom_point).
Removed 3 rows containing missing values (geom_point).
Warning: Removed 2 rows containing missing values (geom_point).
Warning: Removed 3 rows containing missing values (geom_point).

Warning: Removed 1 rows containing missing values (geom_point).
Warning: Removed 2 rows containing missing values (geom_point).
Warning: Removed 1 rows containing missing values (geom_point).
Warning: Removed 2 rows containing missing values (geom_point).

Warning: Removed 1 rows containing missing values (geom_point).
Warning: Removed 1 rows containing missing values (geom_point).

8.1.2 pig skin & skin pad

  # plot skin
for (k in num.var3){
  # get the min/max range of the data set
  range_var <- range(mean.sample[k])  

  skin.mean <- ggplot(data = mean.sample[grep("pig skin", 
               mean.sample[["Contact.material"]]), ], aes_string(x = "Sample", 
               names(mean.sample)[k], colour = "Sample")) + 
               geom_point(size = 4) +
               coord_cartesian(ylim = range_var) +
                 facet_wrap(Sample ~ Contact.material, nrow = 2) + 
               labs(title = "natural contact material") + 
               ylab(names(mean.sample)[j]) + xlab(NULL) +
                 labs(y = gsub("\\.", " ", names(mean.sample)[k])) +
                 scale_colour_futurama() +
                 theme_classic() +
               theme(legend.position = "none") +
               theme(axis.text.x = element_blank(),axis.ticks = element_blank())  
 
   # plot skin pad 
   pad.mean <- ggplot(data = mean.sample[grep("skin pad",
               mean.sample[["Contact.material"]]), ], aes_string(x = "Sample",
               names(mean.sample)[k], colour = "Sample")) + 
                 geom_point(size = 4) +
               coord_cartesian(ylim = range_var) +
                 facet_wrap(Sample ~ Contact.material, nrow = 2) + 
               labs(title = "artificial contact material") + 
               ylab(names(mean.sample)[j]) + xlab(NULL) +
                 labs(y = gsub("\\.", " ", names(mean.sample)[k])) +
                 scale_colour_futurama() +
                 theme_classic() +
               theme(legend.position = "none") +
               theme(axis.text.x = element_blank(),axis.ticks = element_blank())  


    # combine both plots 
  skin.pad.mean <- skin.mean + pad.mean + plot_layout(width = c(3/6, 3/6))  
  print(skin.pad.mean)

  # save to PDF
    file_out <- paste0(file_path_sans_ext(info_in[["file"]]), "_skin_", names(mean.sample)[k],
                ".pdf")
    ggsave(filename = file_out, plot = skin.pad.mean, path = dir_out, device = "pdf", 
           width = 300, height = 170, units = "mm")

} 

Warning: Removed 3 rows containing missing values (geom_point).
Removed 3 rows containing missing values (geom_point).
Removed 3 rows containing missing values (geom_point).
Removed 3 rows containing missing values (geom_point).

Warning: Removed 3 rows containing missing values (geom_point).
Removed 3 rows containing missing values (geom_point).
Removed 3 rows containing missing values (geom_point).
Removed 3 rows containing missing values (geom_point).

Warning: Removed 1 rows containing missing values (geom_point).
Warning: Removed 1 rows containing missing values (geom_point).

8.2 cow scapula & bone plate (with specific range per parameter)

8.2.1 Sq

  # plot cow scapula 
  Sq_bone.mean <- ggplot(data = mean.sample[grep("cow scapula", 
                  mean.sample[["Contact.material"]]), ], aes_string(x = "Sample", 
                  y = "Sq.diff.mean", colour = "Sample")) + 
                  geom_point(size = 4) +
                  coord_cartesian(ylim = c(-2100, 3700)) +
                    facet_wrap(Sample ~ Contact.material, nrow = 2) + 
                  labs(title = "natural contact material") + 
                  ylab("Sq [nm]") + xlab(NULL) +
                    scale_colour_futurama() +
                    theme_classic() +
                  theme(legend.position = "none") +
                  theme(axis.text.x = element_blank(),axis.ticks = element_blank())  
 
   # plot bone plate 
   Sq.plate.mean <- ggplot(data = mean.sample[grep("bone plate",
                    mean.sample[["Contact.material"]]), ], aes_string(x = "Sample", 
                    y = "Sq.diff.mean", colour = "Sample")) + 
                      geom_point(size = 4) +
                    coord_cartesian(ylim = c(-2100, 3700)) +
                      facet_wrap(Sample ~ Contact.material, nrow = 2) + 
                    labs(title = "artificial contact material") + 
                    ylab(NULL) + xlab(NULL) +
                      scale_colour_futurama() +
                      theme_classic() +
                    theme(legend.position = "none") +
                    theme(axis.text.x = element_blank(),axis.ticks = element_blank())  
 
    # combine both plots 
  Sq.bone <- Sq_bone.mean + Sq.plate.mean + plot_layout(width = c(3/6, 3/6))  
  print(Sq.bone)

  # save to PDF
    file_out <- paste0(file_path_sans_ext(info_in[["file"]]), "_Sq.bone.mean", ".pdf")
    ggsave(filename = file_out, plot = Sq.bone, path = dir_out, device = "pdf", 
           width = 300, height = 170, units = "mm")

8.3 cow scapula & bone plate (with specific range per parameter)

8.3.1 Vmc

  # plot cow scapula 
  Vmc_bone.mean <- ggplot(data = mean.sample[grep("cow scapula", 
                   mean.sample[["Contact.material"]]), ], aes_string(x = "Sample", 
                   y = "Vmc.diff.mean", colour = "Sample")) + 
                   geom_point(size = 4) +
                   coord_cartesian(ylim = c(-2, 2.7)) +
                     facet_wrap(Sample ~ Contact.material, nrow = 2) + 
                   labs(title = "natural contact material") + 
                   ylab("Vmc [µm³/µm²]") + xlab(NULL) +
                     scale_colour_futurama() +
                     theme_classic() +
                   theme(legend.position = "none") +
                   theme(axis.text.x = element_blank(),axis.ticks = element_blank())  

   # plot bone plate 
   Vmc.plate.mean <- ggplot(data = mean.sample[grep("bone plate", 
                     mean.sample[["Contact.material"]]), ], aes_string(x = "Sample", 
                     y = "Vmc.diff.mean", colour = "Sample")) + 
                       geom_point(size = 4) +
                     coord_cartesian(ylim = c(-2, 2.7)) +
                       facet_wrap(Sample ~ Contact.material, nrow = 2) + 
                     labs(title = "artificial contact material") + 
                     ylab(NULL) + xlab(NULL) +
                       scale_colour_futurama() +
                       theme_classic() +
                     theme(legend.position = "none") +
                     theme(axis.text.x = element_blank(),axis.ticks = element_blank())  
 
    # combine both plots 
  Vmc.bone <- Vmc_bone.mean + Vmc.plate.mean + plot_layout(width = c(3/6, 3/6))  
  print(Vmc.bone)

  # save to PDF
    file_out <- paste0(file_path_sans_ext(info_in[["file"]]), "_Vmc.bone.mean", ".pdf")
    ggsave(filename = file_out, plot = Vmc.bone, path = dir_out, device = "pdf", 
           width = 300, height = 170, units = "mm")

8.4 cow scapula & bone plate (with specific range per parameter)

8.4.1 HAsfc9

  # plot cow scapula
  HAsfc9_bone.mean <- ggplot(data = mean.sample[grep("cow scapula", 
                      mean.sample[["Contact.material"]]), ], aes_string(x = "Sample", 
                      y = "HAsfc9.diff.mean", colour = "Sample")) + 
                      geom_point(size = 4) +
                      coord_cartesian(ylim = c(-0.5, 10.5)) +
                        facet_wrap(Sample ~ Contact.material, nrow = 2) + 
                      labs(title = "natural contact material") + 
                      ylab("HAsfc9") + xlab(NULL) +
                        scale_colour_futurama() +
                        theme_classic() +
                      theme(legend.position = "none") +
                      theme(axis.text.x = element_blank(),axis.ticks = element_blank())  
                
   # plot bone plate 
   HAsfc9.plate.mean <- ggplot(data = mean.sample[grep("bone plate", 
                        mean.sample[["Contact.material"]]), ], aes_string(x = "Sample", 
                        y = "HAsfc9.diff.mean", colour = "Sample")) + 
                          geom_point(size = 4) +
                        coord_cartesian(ylim = c(-0.5, 10.5)) +
                          facet_wrap(Sample ~ Contact.material, nrow = 2) + 
                        labs(title = "artificial contact material") + 
                        ylab(NULL) + xlab(NULL) +
                          scale_colour_futurama() +
                          theme_classic() +
                        theme(legend.position = "none") +
                        theme(axis.text.x = element_blank(),axis.ticks = element_blank())  

    # combine both plots 
  HAsfc9.bone <- HAsfc9_bone.mean + HAsfc9.plate.mean + plot_layout(width = c(3/6, 3/6))  
  print(HAsfc9.bone)

  # save to PDF
    file_out <- paste0(file_path_sans_ext(info_in[["file"]]), "_HAsfc9.bone.mean", ".pdf")
    ggsave(filename = file_out, plot = HAsfc9.bone, path = dir_out, device = "pdf", 
           width = 300, height = 170, units = "mm")

8.5 cow scapula & bone plate (with specific range per parameter)

8.5.1 epLsar

  # plot cow scapula
  epLsar_bone.mean <- ggplot(data = mean.sample[grep("cow scapula", 
                      mean.sample[["Contact.material"]]), ], aes_string(x = "Sample", 
                      y = "epLsar.diff.mean", colour = "Sample")) + 
                      geom_point(size = 4) +
                      coord_cartesian(ylim = c(-0.002, 0.004)) +
                        facet_wrap(Sample ~ Contact.material, nrow = 2) + 
                      labs(title = "natural contact material") + 
                      ylab("epLsar") + xlab(NULL) +
                        scale_colour_futurama() +
                        theme_classic() +
                      theme(legend.position = "none") +
                      theme(axis.text.x = element_blank(),axis.ticks = element_blank())  
 
   # plot bone plate 
   epLsar.plate.mean <- ggplot(data = mean.sample[grep("bone plate", 
                        mean.sample[["Contact.material"]]), ], aes_string(x = "Sample", 
                        y = "epLsar.diff.mean", colour = "Sample")) + 
                          geom_point(size = 4) +
                        coord_cartesian(ylim = c(-0.002, 0.004)) +
                          facet_wrap(Sample ~ Contact.material, nrow = 2) + 
                        labs(title = "artificial contact material") + 
                        ylab(NULL) + xlab(NULL) +
                          scale_colour_futurama() +
                          theme_classic() +
                        theme(legend.position = "none") +
                        theme(axis.text.x = element_blank(),axis.ticks = element_blank())  

    # combine both plots 
  epLsar.bone <- epLsar_bone.mean + epLsar.plate.mean + plot_layout(width = c(3/6, 3/6))  
  print(epLsar.bone)

  # save to PDF
    file_out <- paste0(file_path_sans_ext(info_in[["file"]]), "_epLsar.bone.mean", ".pdf")
    ggsave(filename = file_out, plot = epLsar.bone, path = dir_out, device = "pdf", 
           width = 300, height = 170, units = "mm")

8.6 cow scapula & bone plate (with specific range per parameter)

8.6.1 Asfc

  # plot cow scapula
  Asfc_bone.mean <- ggplot(data = mean.sample[grep("cow scapula", 
                    mean.sample[["Contact.material"]]), ], aes_string(x = "Sample", 
                    y = "Asfc.diff.mean", colour = "Sample")) + 
                    geom_point(size = 4) +
                    coord_cartesian(ylim = c(-16, 26)) +
                      facet_wrap(Sample ~ Contact.material, nrow = 2) + 
                    labs(title = "natural contact material") + 
                    ylab("Asfc") + xlab(NULL) +
                      scale_colour_futurama() +
                      theme_classic() +
                    theme(legend.position = "none") +
                    theme(axis.text.x = element_blank(),axis.ticks = element_blank())  
 
   # plot bone plate 
   Asfc.plate.mean <- ggplot(data = mean.sample[grep("bone plate", 
                      mean.sample[["Contact.material"]]), ], aes_string(x = "Sample", 
                      y = "Asfc.diff.mean", colour = "Sample")) + 
                        geom_point(size = 4) +
                      coord_cartesian(ylim = c(-16, 26)) +
                        facet_wrap(Sample ~ Contact.material, nrow = 2) + 
                      labs(title = "artificial contact material") + 
                      ylab(NULL) + xlab(NULL) +
                        scale_colour_futurama() +
                        theme_classic() +
                      theme(legend.position = "none") +
                      theme(axis.text.x = element_blank(),axis.ticks = element_blank())  

    # combine both plots 
  Asfc.bone <- Asfc_bone.mean + Asfc.plate.mean + plot_layout(width = c(3/6, 3/6))  
  print(Asfc.bone)

  # save to PDF
    file_out <- paste0(file_path_sans_ext(info_in[["file"]]), "_Asfc.bone.mean", ".pdf")
    ggsave(filename = file_out, plot = Asfc.bone, path = dir_out, device = "pdf", 
           width = 300, height = 170, units = "mm")

8.7 pig skin & skin pad (with specific range per parameter)

8.7.1 Sq

  # plot pig skin 
  Sq_skin.mean <- ggplot(data = mean.sample[grep("pig skin", 
                  mean.sample[["Contact.material"]]), ], aes_string(x = "Sample",  
                  y = "Sq.diff.mean", colour = "Sample")) + 
                  geom_point(size = 4) +
                  coord_cartesian(ylim = c(-2100, 3700)) +
                    facet_wrap(Sample ~ Contact.material, nrow = 2) + 
                  labs(title = "natural contact material") + 
                  ylab("Sq [nm]") + xlab(NULL) +
                    scale_colour_futurama() +
                    theme_classic() +
                  theme(legend.position = "none") +
                  theme(axis.text.x = element_blank(),axis.ticks = element_blank())  
                
   # plot skin pad 
   Sq.pad.mean <- ggplot(data = mean.sample[grep("skin pad", 
                  mean.sample[["Contact.material"]]), ], aes_string(x = "Sample", 
                  y = "Sq.diff.mean", colour = "Sample")) + 
                    geom_point(size = 4) +
                  coord_cartesian(ylim = c(-2100, 3700)) +
                    facet_wrap(Sample ~ Contact.material, nrow = 2) + 
                  labs(title = "artificial contact material") + 
                  ylab(NULL) + xlab(NULL) +
                    scale_colour_futurama() +
                    theme_classic() +
                  theme(legend.position = "none") +
                  theme(axis.text.x = element_blank(),axis.ticks = element_blank())  

    # combine both plots 
  Sq.skin <- Sq_skin.mean + Sq.pad.mean + plot_layout(width = c(3/6, 3/6))  
  print(Sq.skin)

  # save to PDF
    file_out <- paste0(file_path_sans_ext(info_in[["file"]]), "_Sq.skin.mean", ".pdf")
    ggsave(filename = file_out, plot = Sq.skin, path = dir_out, device = "pdf", 
           width = 300, height = 170, units = "mm")

8.8 pig skin & skin pad (with specific range per parameter)

8.8.1 Vmc

  # plot pig skin 
  Vmc_skin.mean <- ggplot(data = mean.sample[grep("pig skin", 
                   mean.sample[["Contact.material"]]), ], aes_string(x = "Sample", 
                   y = "Vmc.diff.mean", colour = "Sample")) + 
                   geom_point(size = 4) +
                   coord_cartesian(ylim = c(-2, 2.7)) +
                     facet_wrap(Sample ~ Contact.material, nrow = 2) + 
                   labs(title = "natural contact material") + 
                   ylab("Vmc [µm³/µm²]") + xlab(NULL) +
                     scale_colour_futurama() +
                     theme_classic() +
                   theme(legend.position = "none") +
                   theme(axis.text.x = element_blank(),axis.ticks = element_blank())  
    
   # plot skin pad 
   Vmc.pad.mean <- ggplot(data = mean.sample[grep("skin pad", 
                   mean.sample[["Contact.material"]]), ], aes_string(x = "Sample", 
                   y = "Vmc.diff.mean", colour = "Sample")) + 
                     geom_point(size = 4) +
                   coord_cartesian(ylim = c(-2, 2.7)) +
                     facet_wrap(Sample ~ Contact.material, nrow = 2) + 
                   labs(title = "artificial contact material") + 
                   ylab(NULL) + xlab(NULL) +
                     scale_colour_futurama() +
                     theme_classic() +
                   theme(legend.position = "none") +
                   theme(axis.text.x = element_blank(),axis.ticks = element_blank())  

    # combine both plots 
  Vmc.skin <- Vmc_skin.mean + Vmc.pad.mean + plot_layout(width = c(3/6, 3/6))  
  print(Vmc.skin)

  # save to PDF
    file_out <- paste0(file_path_sans_ext(info_in[["file"]]), "_Vmc.skin.mean", ".pdf")
    ggsave(filename = file_out, plot = Vmc.skin, path = dir_out, device = "pdf", 
           width = 300, height = 170, units = "mm")

8.9 pig skin & skin pad (with specific range per parameter)

8.9.1 HAsfc9

  # plot pig skin 
  HAsfc9_skin.mean <- ggplot(data = mean.sample[grep("pig skin", 
                      mean.sample[["Contact.material"]]), ], aes_string(x = "Sample", 
                      y = "HAsfc9.diff.mean", colour = "Sample")) + 
                      geom_point(size = 4) +
                      coord_cartesian(ylim = c(-0.5, 10.5)) +
                        facet_wrap(Sample ~ Contact.material, nrow = 2) + 
                      labs(title = "natural contact material") + 
                      ylab("HAsfc9") + xlab(NULL) +
                        scale_colour_futurama() +
                        theme_classic() +
                      theme(legend.position = "none") +
                      theme(axis.text.x = element_blank(),axis.ticks = element_blank())  
                
   # plot skin pad 
   HAsfc9.pad.mean <- ggplot(data = mean.sample[grep("skin pad", 
                      mean.sample[["Contact.material"]]), ], aes_string(x = "Sample", 
                      y = "HAsfc9.diff.mean", colour = "Sample")) + 
                        geom_point(size = 4) +
                      coord_cartesian(ylim = c(-0.5, 10.5)) +
                        facet_wrap(Sample ~ Contact.material, nrow = 2) + 
                      labs(title = "artificial contact material") + 
                      ylab(NULL) + xlab(NULL) +
                        scale_colour_futurama() +
                        theme_classic() +
                      theme(legend.position = "none") +
                      theme(axis.text.x = element_blank(),axis.ticks = element_blank())  

    # combine both plots 
  HAsfc9.skin <- HAsfc9_skin.mean + HAsfc9.pad.mean + plot_layout(width = c(3/6, 3/6))   
  print(HAsfc9.skin)

  # save to PDF
    file_out <- paste0(file_path_sans_ext(info_in[["file"]]), "_HAsfc9.skin.mean", ".pdf")
    ggsave(filename = file_out, plot = HAsfc9.skin, path = dir_out, device = "pdf", 
           width = 300, height = 170, units = "mm")

8.10 pig skin & skin pad (with specific range per parameter)

8.10.1 epLsar

  # plot pig skin 
  epLsar_skin.mean <- ggplot(data = mean.sample[grep("pig skin",  
                      mean.sample[["Contact.material"]]), ], aes_string(x = "Sample", 
                      y = "epLsar.diff.mean", colour = "Sample")) + 
                      geom_point(size = 4) +
                      coord_cartesian(ylim = c(-0.002, 0.004)) +
                        facet_wrap(Sample ~ Contact.material, nrow = 2) + 
                      labs(title = "natural contact material") + 
                      ylab("epLsar") + xlab(NULL) +
                        scale_colour_futurama() +
                        theme_classic() +
                      theme(legend.position = "none") +
                      theme(axis.text.x = element_blank(),axis.ticks = element_blank())  
                
   # plot skin pad 
   epLsar.pad.mean <- ggplot(data = mean.sample[grep("skin pad", 
                      mean.sample[["Contact.material"]]), ], aes_string(x = "Sample", 
                      y = "epLsar.diff.mean", colour = "Sample")) + 
                        geom_point(size = 4) +
                      coord_cartesian(ylim = c(-0.002, 0.004)) +
                        facet_wrap(Sample ~ Contact.material, nrow = 2) + 
                      labs(title = "artificial contact material") + 
                      ylab(NULL) + xlab(NULL) +
                        scale_colour_futurama() +
                        theme_classic() +
                      theme(legend.position = "none") +
                      theme(axis.text.x = element_blank(),axis.ticks = element_blank())  

    # combine both plots 
  epLsar.skin <- epLsar_skin.mean + epLsar.pad.mean + plot_layout(width = c(3/6, 3/6))  
  print(epLsar.skin)

  # save to PDF
    file_out <- paste0(file_path_sans_ext(info_in[["file"]]), "_epLsar.skin.mean", ".pdf")
    ggsave(filename = file_out, plot = epLsar.skin, path = dir_out, device = "pdf", 
            width = 300, height = 170, units = "mm")

8.11 pig skin & skin pad (with specific range per parameter)

8.11.1 Asfc

  # plot pig skin 
  Asfc_skin.mean <- ggplot(data = mean.sample[grep("pig skin", 
                    mean.sample[["Contact.material"]]), ], aes_string(x = "Sample", 
                    y = "Asfc.diff.mean", colour = "Sample")) + 
                    geom_point(size = 4) +
                    coord_cartesian(ylim = c(-16, 26)) +
                      facet_wrap(Sample ~ Contact.material, nrow = 2) + 
                    labs(title = "natural contact material") + 
                    ylab("Asfc") + xlab(NULL) +
                      scale_colour_futurama() +
                      theme_classic() +
                    theme(legend.position = "none") +
                    theme(axis.text.x = element_blank(),axis.ticks = element_blank())  
                
   # plot skin pad 
   Asfc.pad.mean <- ggplot(data = mean.sample[grep("skin pad", 
                    mean.sample[["Contact.material"]]), ], aes_string(x = "Sample", 
                    y = "Asfc.diff.mean", colour = "Sample")) + 
                      geom_point(size = 4) +
                    coord_cartesian(ylim = c(-16, 26)) +
                      facet_wrap(Sample ~ Contact.material, nrow = 2) + 
                    labs(title = "artificial contact material") + 
                    ylab(NULL) + xlab(NULL) +
                      scale_colour_futurama() +
                      theme_classic() +
                    theme(legend.position = "none") +
                    theme(axis.text.x = element_blank(),axis.ticks = element_blank())  

    # combine both plots 
  Asfc.skin <- Asfc_skin.mean + Asfc.pad.mean + plot_layout(width = c(3/6, 3/6))  
  print(Asfc.skin)

  # save to PDF
    file_out <- paste0(file_path_sans_ext(info_in[["file"]]), "_Asfc.skin.mean", ".pdf")
    ggsave(filename = file_out, plot = Asfc.skin, path = dir_out, device = "pdf", 
           width = 300, height = 170, units = "mm")

9 Principal component analysis

9.1 PCA contact material - raw material seperated

9.1.1 flint

# filter the data for flint only
flint <- filter(imp_data, Raw.material == "flint") 

# remove rows with na values 
data_pca.flint <- na.omit(flint)
# use for the PCA only selected variables: Sq, SSK, Vmc, Isotropy, Mean density of furrows,
# Asfc, HAsfc9 
flint.pca <- prcomp(data_pca.flint[ , c(9:10, 27, 32:33, 40, 42)], scale. = TRUE) 
# convert the data into factor 
data_pca.flint[["Contact.material"]] <- factor(data_pca.flint[["Contact.material"]])

# using ggfortify
PCA.flint <- autoplot(flint.pca, data = data_pca.flint, colour = "Contact.material", 
             size = 2,
             loadings = TRUE, loadings.colour = "black", loadings.label = TRUE,
             loadings.label.colour = "black", 
             loadings.label.size  = 4, loadings.label.repel = TRUE,  
             frame.alpha = 0) + 
             theme_classic() +
             scale_colour_futurama() 
         
print(PCA.flint)

# save the plot
file_out <- paste0(file_path_sans_ext(info_in[["file"]]), "_PCA.flint_contact", ".pdf")
ggsave(filename = file_out, plot = PCA.flint, path = dir_out, device = "pdf")

9.1.2 lydite

# filter the data for lydite only
lydite <- filter(imp_data, Raw.material == "lydite") 

# remove rows with na values 
data_pca.lydite <- na.omit(lydite)
# use for the PCA only selected variables: Sq, SSK, Vmc, Isotropy, Mean density of furrows,
# Asfc, HAsfc9 
lydite.pca <- prcomp(data_pca.lydite[ , c(9:10, 27, 32:33, 40, 42)], scale. = TRUE) 
# convert the data into factor 
data_pca.lydite[["Contact.material"]] <- factor(data_pca.lydite[["Contact.material"]])

# using ggfortify
PCA.lydite <- autoplot(lydite.pca, data = data_pca.lydite, colour = "Contact.material", 
              size = 2,
              loadings = TRUE, loadings.colour = "black", loadings.label = TRUE,
              loadings.label.colour = "black", 
              loadings.label.size  = 4, loadings.label.repel = TRUE,  
              frame.alpha = 0) + 
              theme_classic() +
              scale_colour_futurama() 
         
print(PCA.lydite)

# save the plot
file_out <- paste0(file_path_sans_ext(info_in[["file"]]), "_PCA.lydite_contact", ".pdf")
ggsave(filename = file_out, plot = PCA.lydite, path = dir_out, device = "pdf")

9.2 PCA before and after - raw material seperated

9.2.1 flint

# using ggfortify
PCA2_flint<- autoplot(flint.pca, data = data_pca.flint, colour = "Cycle", size = 2,
             loadings = TRUE, loadings.colour = "black", loadings.label = TRUE, 
             loadings.label.colour = "black", 
             loadings.label.size  = 4, loadings.label.repel = TRUE,  
             frame = TRUE, frame.type = "convex", frame.colour = "Cycle", frame.alpha = 0) +
             theme_classic() +
             scale_colour_manual(values = custom.col6$col)
Error in is_missing(values): object 'custom.col6' not found
print(PCA2_flint)
Error in print(PCA2_flint): object 'PCA2_flint' not found
# save the plot
file_out <- paste0(file_path_sans_ext(info_in[["file"]]), "_PCA.flint_cycle", ".pdf")
ggsave(filename = file_out, plot = PCA2_flint, path = dir_out, device = "pdf")
Error in plot_theme(plot): object 'PCA2_flint' not found

9.2.2 lydite

# using ggfortify
PCA2_lydite<- autoplot(lydite.pca, data = data_pca.lydite, colour = "Cycle", size = 2,
              loadings = TRUE, loadings.colour = "black", loadings.label = TRUE, 
              loadings.label.colour = "black", 
              loadings.label.size  = 4, loadings.label.repel = TRUE,  
              frame = TRUE, frame.type = "convex", frame.colour = "Cycle", 
              frame.alpha = 0) +
              theme_classic() +
              scale_colour_manual(values = custom.col6$col)
Error in is_missing(values): object 'custom.col6' not found
print(PCA2_lydite)
Error in print(PCA2_lydite): object 'PCA2_lydite' not found
# save the plot
file_out <- paste0(file_path_sans_ext(info_in[["file"]]), "_PCA.lydite_cycle", ".pdf")
ggsave(filename = file_out, plot = PCA2_lydite, path = dir_out, device = "pdf")
Error in plot_theme(plot): object 'PCA2_lydite' not found

10 sessionInfo() and RStudio version

sessionInfo()
R version 4.1.0 (2021-05-18)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 10 x64 (build 19043)

Matrix products: default

locale:
[1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252   
[3] LC_MONETARY=German_Germany.1252 LC_NUMERIC=C                   
[5] LC_TIME=German_Germany.1252    

attached base packages:
[1] tools     stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] ggfortify_0.4.14  ggsci_2.9         patchwork_1.1.1   forcats_0.5.1    
 [5] stringr_1.4.0     dplyr_1.0.9       purrr_0.3.4       readr_2.1.2      
 [9] tidyr_1.2.0       tibble_3.1.6      tidyverse_1.3.1   doBy_4.6.13      
[13] ggplot2_3.3.6     R.utils_2.11.0    R.oo_1.24.0       R.methodsS3_1.8.1

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.8.3         ggrepel_0.9.1        lubridate_1.8.0     
 [4] lattice_0.20-44      assertthat_0.2.1     digest_0.6.29       
 [7] utf8_1.2.2           R6_2.5.1             cellranger_1.1.0    
[10] backports_1.4.1      reprex_2.0.1         evaluate_0.15       
[13] highr_0.9            httr_1.4.2           pillar_1.7.0        
[16] rlang_1.0.2          readxl_1.4.0         rstudioapi_0.13     
[19] jquerylib_0.1.4      Matrix_1.3-3         rmarkdown_2.14      
[22] labeling_0.4.2       munsell_0.5.0        broom_0.8.0         
[25] compiler_4.1.0       Deriv_4.1.3          modelr_0.1.8        
[28] xfun_0.30            pkgconfig_2.0.3      microbenchmark_1.4.9
[31] htmltools_0.5.2      tidyselect_1.1.2     gridExtra_2.3       
[34] fansi_1.0.3          crayon_1.5.1         tzdb_0.3.0          
[37] dbplyr_2.1.1         withr_2.5.0          MASS_7.3-54         
[40] grid_4.1.0           jsonlite_1.8.0       gtable_0.3.0        
[43] lifecycle_1.0.1      DBI_1.1.2            magrittr_2.0.3      
[46] scales_1.2.0         cli_3.3.0            stringi_1.7.6       
[49] farver_2.1.0         fs_1.5.2             xml2_1.3.3          
[52] bslib_0.3.1          ellipsis_0.3.2       generics_0.1.2      
[55] vctrs_0.4.1          glue_1.6.2           hms_1.1.1           
[58] fastmap_1.1.0        yaml_2.3.5           colorspace_2.0-3    
[61] rvest_1.0.2          knitr_1.39           haven_2.5.0         
[64] sass_0.4.1          

RStudio version 1.4.1717.


11 Cite R packages used

for (i in pack_to_load) print(citation(i), bibtex = FALSE)

To cite package 'R.utils' in publications use:

  Henrik Bengtsson (2021). R.utils: Various Programming Utilities. R
  package version 2.11.0. https://CRAN.R-project.org/package=R.utils


The 'tools' package is part of R.  To cite R in publications use:

  R Core Team (2021). R: A language and environment for statistical
  computing. R Foundation for Statistical Computing, Vienna, Austria.
  URL https://www.R-project.org/.

We have invested a lot of time and effort in creating R, please cite it
when using it for data analysis. See also 'citation("pkgname")' for
citing R packages.


To cite ggplot2 in publications, please use:

  H. Wickham. ggplot2: Elegant Graphics for Data Analysis.
  Springer-Verlag New York, 2016.


To cite package 'doBy' in publications use:

  Søren Højsgaard and Ulrich Halekoh (2022). doBy: Groupwise
  Statistics, LSmeans, Linear Estimates, Utilities. R package version
  4.6.13. https://CRAN.R-project.org/package=doBy

ATTENTION: This citation information has been auto-generated from the
package DESCRIPTION file and may need manual editing, see
'help("citation")'.


  Wickham et al., (2019). Welcome to the tidyverse. Journal of Open
  Source Software, 4(43), 1686, https://doi.org/10.21105/joss.01686


To cite package 'patchwork' in publications use:

  Thomas Lin Pedersen (2020). patchwork: The Composer of Plots. R
  package version 1.1.1. https://CRAN.R-project.org/package=patchwork


To cite package 'ggsci' in publications use:

  Nan Xiao (2018). ggsci: Scientific Journal and Sci-Fi Themed Color
  Palettes for 'ggplot2'. R package version 2.9.
  https://CRAN.R-project.org/package=ggsci


To cite ggfortify in publications, please use:

  Yuan Tang, Masaaki Horikoshi, and Wenxuan Li. "ggfortify: Unified
  Interface to Visualize Statistical Result of Popular R Packages." The
  R Journal 8.2 (2016): 478-489.

  Masaaki Horikoshi and Yuan Tang (2016). ggfortify: Data Visualization
  Tools for Statistical Analysis Results.
  https://CRAN.R-project.org/package=ggfortify

END OF SCRIPT